Force distribution in granular media studied by an energy method based on statistical 

mechanics 
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In the present letter a method to find a proper expression for the force distribution inside a 
granular sample in static equilibrium is proposed. The method is based in statistical mechanics and 
the force distribution is obtained by studying how the potential elastic energy is divided among 
the different contacts between grains. It is found with DEM simulations with spheres that the 
elastic potential energy distribution follows a Bose Einstein law from which the force distribution is 
deduced. The present letter open a way in which granular materials can be studied with the tools 
provided by statistical mechanics. 



The connection between the macroscopical behavior of 
granular materials and their microscopic interactions is 
still an open question with many unsolved problems. The 
main difficulty is that the methods proposed in statisti- 
cal mechanics, the main branch of physics studying the 
micro-macro connection, cannot be easily used directly in 
granular materials due to their great energy dissipation. 
Considering this problems, Edwards and BumenfeldfH 
proposed a partition function based on the volume occu- 
pied by the grains and no on their individual energy. This 
certainly avoids the energy dissipation problems and also 
puts the volume of the ensemble as the main state vari- 
able of the granular ensemble. However the computation 
of the number of states associated with a given volume 
can be difficult because of many factors like for example 
the complex shape of particles. Since each particle has 
six degrees of freedom, the integral on the configuration 
space can be quite complex. Moreover the volumes of 
the particles interact with each other since the total vol- 
ume must remain constant (as the energy in an ideal gas) 
which add even more complexity to the calculation of the 
configuration integrals. 

With DEM based simulations the force statistic have 
been studied and it has been found by Radjai, with Con- 
tact Dynamics simulations, that the distribution of the 
values of normal forces behaves like a power law be- 
fore the mean value and like a exponential decay after- 
wards 0], however there is currently no justification for 
this particular behavior yet. Liu et. [3] provided a new 
force distribution based on a mean field approximation 
called the g- model, but still it fails to properly represent 
small forces presented in the ensemble as discussed by 
Luding et. al\^. 

In the present letter the author shows some results 
of force statistics obtained by DEM simulations with 
spheres in 3D. The force distribution is deduced by find- 
ing the distribution of potential energy in the contacts. 
The contacts are taken then as the new particles from 
which the ensemble is built. This particles, called by 
the author Elastons, represent the elastic interaction 
between two grains and as will be shown they greatly 
simplify the statistical mechanical study of granular ma- 



terials. 

The proposed model Although highly dissipative, 
granular materials follow the first law of thermodynamics 
as any other physical entity. 



dU = 5Q- SW, 



(1) 



with U the change in internal energy, Q the heat dissi- 
pated by the system and W the mechanical work done 
by the system. In granular materials this should be the 
template for every constitute model since the internal en- 
ergy represents the potential energy stored by the elastic 
properties of the material, the heat represent the energy 
dissipated by frictional and viscous forces and hence the 
plastic deformation, and the work done relates the ex- 
ternal stresses and strains. However given the difficulty 
in finding the functional form of each of these terms, the 
first law of thermodynamics is rarely the foundation of 
constitutive laws. In the present letter the author will 
focus his attention to the term related with the elastic 
domain: The internal energy U that is not dissipated by 
viscous or frictional forces. 

Once the kinetic energy of the grains is dissipated by 
internal viscosity and friction, the only mechanical energy 
remaining is the elastic potential energy stored in the 
contact points. This contact points are discrete and the 
total internal energy is just the addition of the individual 
potential energies at each contact point. These Elastons 
distribute the internal energy U in the granular ensemble 
among themselves. Now in a granular ensemble of N 
spherical grains there can be as much as N{N — l)/2 
possible Elastons but most of them are not relevant since 
a spherical grain is only in contact with just a few other 
grains in general. However Elastons with non-zero energy 
are constantly being created and destroyed as the sample 
is subjected to external stresses representing the change 
in structure. 

In the initial state and in absence of external forces, the 
grains are just in contact with each other but there is no 
repulsion force. We called this state the ground state and 
all the Elastons have the same energy equal to zero. Now 
if stresses are applied to the sample some contacts start 
to accumulate potential energy and hence some Elastons 
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jump from the ground state to higher energy levels. At 
the end of a compression test, most Elastons will have 
the lowest energy possible while some others will have 
higher energies. Considering that Elastons are in princi- 
ple statistically indistinguishably, is tempting to assign a 
Bose-Einsten distribution form for their potential energy 
E, 



P{E) 



9{E) 



cxp(/3(i? - /i)) - r 



(2) 



with /3 related to the mean energy of the ensemble, /i a 
parameter that in the case of gases represent the chemical 
potential and g{E) the density of states. Now, is danger- 
ous to immediately assign the same physical meaning of 
these quantities for the granular material as for the par- 
ticle gases studied in statistical mechanics. For example 
the factor /3 is related to the temperature of the gas as 
measured by a thermometer. But we cannot measure the 
P factor in the same way as with the gases since in a gas 
the temperature is measured by radiation or convection, 
i. e. by the movement of particles. The Elastons do 
not move, their energy is entirely potential energy and 
hence the thermometer will not give us an accurate mea- 
surement of their granular temperature. At most, we can 
establish (3 and /i as parameters to be found by fitting 
the above relation with experimental data. 

But how to measure the potential energy of an Elaston? 
Clearly is very difficult to do so in an experimental setup. 
However the DEM method introduced by Cundall Q for 
the study of granular materials allows us to measure these 
microscopic quantities with ease. All that is needed is a 
contact law relating the potential energy and the force 
exerted by the grains. 

In order to simplify the mathematical deductions the 
author has chosen the linear spring law that has been 
used recently for the simulation of the True Triaxial Test 
by Donze et. al.^ instead of the more popular Hertz 
law . In the linear spring law the spheres are assumed 
to overlap with each other and the elastic force is pro- 
portional to the overlapping length [§1 . The form of this 
force depends on a stiffness constant Kn as. 



F = KJ, 



(3) 



where 5 — Sh with n the normal vector joining the two 
sphere's centers and d is the overlapping length. 

With these simple model for the contact force, the en- 
ergy of the corresponding Elaston can be found. 



E = ^KJ^ 



(4) 



In order to reach the static equilibrium state the whole 
kinetic energy has to be dissipated by viscous forces as 
used in previous studies Once the static equilibrium 
state is reached, the internal energy is distributed among 
the Elastons as mentioned above. To apply Eq. [2]to the 



distribution of potential energy, the density of states for 
a given energy must be calculated. Elastons have only 
three degrees of freedom given by the components of the 
overlapping vector S. The number of possible states 
for an Elaston with overlapping length smaller than 6 is 
then proportional to the volume integral in the configu- 
ration space. 



dSxdSydS: 



= J AnS^dS = (5) 



And hence the number of states with an energy smaller 
than E is then by Eq. lU 



, , Att f2E 



(6) 



Therefore the density of states can easily be found by 
differentiation of the previous function with respect to 
E, 



, , d^(E) 47r / 2 \ ^ ^ 



dE " 2 \KnJ " ' '■''^ 
and accordingly Eq. [5]is rewritten in the present case as: 

e'2 



P{E) oc 



exp(/3(i?-M))-l- 



(8) 



Equations [S] and [3] can be used to obtain the distri- 
bution of forces P{F) in the static equilibrium state of 
the granular ensemble. In order to do this, the reader 
must remember that for two random variables x and y 
related by a function y = <^{x) their distributions follow 



P{y) = P{x) 



and consequently. 



P{F) oc 



exp 



/i 



(9) 



Simulations In order to check the distribution laws 
several True Triaxial Test (TTT) simulation have been 
carried out. In the TTT setup the spheres are enclosed 
by a set of six walls. Forces are applied to the walls 
in order to obtain a constant hydrostatic pressure over 
the sample. The simulation runs until the kinetic en- 
ergy is dissipated. Once the static equilibrium state is 
reached the potential energy at each contact is measured 
and this is the energy associated to the Elaston. In the 
present letter some randomness is introduce by setting 
the spheres initially at a Hexagonal Close Packing and 
then randomly deleting some of them until the desired 
number of spheres is met. A total of 8 different hydro- 
static pressures were considered with 8 random samples 
per case. 
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Results For a given pressure the results for the force 
distribution at the static equihbrium test are shown in 
Fig. [T] As can be seen the fitting is quite good with 
a i?2 factor of 0.999 with values /? = 5.13 ± 0.12J^i 
and /Li = —0.0056 ± 0.0009J. An important conclusion 
form the fitting process is that n must be smaller than 
zero, otherwise there can be force values with negative 
population. In the simulation data the behavior observed 
by Radjai [1] is also present as the force distribution can 
be approximated by a power law below certain value and 
by a exponential law after certain threshold. The model 
proposed in this letter does not consider this especial 
cases and covers the whole range of forces. 
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FIG. 1: Normalized force histogram (circles) for a given pres- 
sure in a) Linear-Linear scale, b) Linear Log scale and c) Log 
Log scale. The results are compared with the model of Eq. [9] 
(solid line). 

As the energy is the base of the proposed model, it was 
also fitted with Eq. [51 The results are shown in Fig. [2] 
where the same values for /3 and fi from the previous 
fitting have been used leaving only the proportionality 
constant as the only fitting parameter. Again the fitting 
has a good value of 0.999. As expected the majority 
of Elusions lie in levels of small energies. It was also 
compared with the Maxwell Boltzmann statistics widely 
used for the macroscopic scale in which the population 
of energy states is given by: 



Clearly by observing the DEM data, the BE model is a 
better fit than the classical MB model, confirming again 
that the Elastons behave statistically like Bosons. 
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FIG. 2: Normalized histogram of potential energy occupation 
(circles) for a given pressure in a) Linear-Linear scale, b) Log 
Log scale. The results are compared with the BE model of 
Eq. [8]and also with the Maxwell Boltzmann (MB) model. 

Next the dependence of the potential energy distribu- 
tion on the applied hydrostatic pressure is observed. In 
Fig. [3] the results for three different pressures are shown. 
The process of applying confining pressure to the Triaxial 
Test gives energy to the whole system and hence Elas- 
tons leave the ground state and occupy levels with higher 
energies. As expected the minimum pressure considered 
guaranties a low dispersion of Elastons over the energy 
spectrum, and the highest pressure ensures the largest 
dispersion. 

In a gas of particles, the dispersion of particles across 
the energy spectrum is usually related to the tempera- 
ture. In the case of the present letter, the analogy should 
be done with the /? parameter of Eqs. [5]and[ni Hence this 
parameter must be checked against the pressure applied 
to the confinement. This is done by plotting the inverse 
of (3 against the work done on the sample expressed as 
W — P/SV where /SV is the difference in volume be- 
tween the first pressure considered p = 1 to the current 
pressure. In Fig. 



4(a) it is shown how the inverse of 



P{E) oc g{E) exp {~/3E) (x E^ exp (—/?£') (10) /3 grows with the external work done to the system by 
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FIG. 3: Normalized energy histogram for three different pres- 
sures (expressed in Pa) in Log-Linear scale for detailed obser- 
vation. 

the pressure. It is not a linear relation, and its functional 
form is difficult to establish at the moment since the work 
done is divided in the potential elastic energy among the 
Elastons and in the dissipated energy that is not being 
considered here. 
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FIG. 4: Parameter /3 versus a) the work done to the system 
PAV and b) The mean energy of the Elastons. 

The mean energy per Elaston is in contrast propor- 
tional to the inverse of (3 (Fig. [4(b)] ) as the mean en- 



ergy of particles is proportional to the temperature in a 
gas. This is indeed expected an open a possible way to 
measure the mean potential energy by a special granular 
temperature. 

Conclusions In the present letter a method to study 
the statistics of force chains inside granular materials is 
proposed, the method is based in energy considerations. 
The question of how the potential elastic energy is dis- 
tributed among the contacts is solved by a Bose Einstein 
(BE) law rather than a Maxwell Boltzmann (MB) one 
and, thanks to this, a force distribution is obtained and 
checked with DEM data. This implies that the contacts 
can be represented by particles that the author call Elas- 
tons and the tools provided by statistical mechanics for 
the study of Bosons gases can be use for the granular gas. 
The Elastons are far better suited for study since they 
have only three degrees of freedom in contrast with the 
six degrees that a grain in 3D space has. Moreover this 
particles do not interact with each other, which greatly 
simplifies the calculation of statistical properties. Future 
research in this area should focus in how the parame- 
ters of the proposed model area affected by the external 
stresses of the granular material, and that will help to the 
formulation of constitutive laws based on first principles. 
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